arXivrl 506.05473V 1 [astro-ph.HE] 17Jun2015 


Black hole-neutron star binary merger: Dependence on black hole spin orientation 

and equations of state 

Kyohei Kawaguchi/ Koutarou Kyutoku/ Hiroyuki Nakano/’^ 

Hirotada Okawa/’^ Masam Shibata/ and Keisuke Taniguchi^’^ 

^ Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan 
‘^Interdisciplinary Theoretical Science (iTHES) Research Group, RIKEN, Wako, Saitama 351-0198, Japan 
^Department of Physics, Kyoto University, Kyoto 606-8502, Japan 
^Center for Computational Relativity and Gravitation, 
and School of Mathematical Sciences, Rochester Institute of Technology, 

85 Lomh Memorial Drive, Rochester, New York 14623 
^Advanced Research Institute for Science and Engineering, 

Waseda University, 3-4-1 Okuho, Shinjuku, Tokyo 169-8555, Japan 
^Department of Physics, University of the Ryukyus, 

1 Senbaru, Nishihara, Okinawa 903-0213, Japan 

^Graduate School of Arts and Sciences, University of Tokyo, Komaha, Meguro, Tokyo 153-8902, Japan 

(Dated: June 19, 2015) 

We systematically performed numerical-relativity simulations for black hole (BH) - neutron star 
(NS) binary mergers with a variety of the BH spin orientation and nuclear-theory-based equations of 
state (EOS) of the NS. The initial misalignment angles of the BH spin measured from the direction 
of the orbital angular momentum are chosen in the range of hiit,o ~ 30°-90°. We employed four 
models of nuclear-theory-based zero-temperature EOS for the NS with which the compactness of 
the NS is in the range of C = Mns/J?ns = 0.138-0.180, where Mns and i^NS are the mass and 
the radius of the NS, respectively. The mass ratio of the BH to the NS, Q = Mbh/Mns, and the 
dimensionless spin parameter of the BH, y, are chosen to be Q = 5 and y = 0.75, together with 
TTns = 1.35Mq so that the BH spin misalignment has a significant effect on tidal disruption of the 
NS. We obtain the following results: (i) The inclination angle of hiit,o < 70° and hiit,o < 50° are 
required for the formation of a remnant disk with its mass larger than O.IM© for the case C = 0.140 
and C — 0.160, respectively, while the disk mass is always smaller than O.IM© for C > 0.175. The 
ejecta with its mass larger than O.OIM© is obtained for Ztiit,o < 85° with C = 0.140, for Ztiit,o < 65° 
with C — 0.160, and for hiit,o < 30° with C — 0.175. (ii) The rotational axis of the dense part of the 
remnant disk with its rest-mass density larger than 10® g/cm^ is approximately aligned with the 
remnant BH spin for ztiit,o ~ 30°. On the other hand, the disk axis is misaligned initially with ~ 30° 
for hiit,o ~ 60°, and the alignment with the remnant BH spin is achieved at ~ 50-60 ms after the 
onset of merger. The accretion time scale of the remnant disk is typically ~ 100 ms and depends 
only weakly on the misalignment angle and the EOS. (iii) The ejecta velocity is typically ~ 0.2-0.3c 
and depends only weakly on the misalignment angle and the EOS of the NS, while the morphology 
of the ejecta depends on its mass, (iv) The gravitational-wave spectra contains the information of 
the NS compactness in the cutoff frequency for hiit,o ^ 60°. 

PACS numbers: 04.25.D-, 04.30.-w, 04.40.Dg 

I. INTRODUCTION 

The ground-based gravitational-wave detectors, such 
as Advanced LIGO [T], Advanced VIRGO [2], and KA- 
GRA [3], will detect signals in the next decade. Black 
hole (BH) - neutron star (NS) binary mergers are one of 
the most promising gravitational-wave sources for these 
detectors [4]. Since NSs and BHs are compact objects 
with strong gravity, detection of gravitational waves from 
compact binary mergers, including BH-NS mergers, can 
be a touchstone of the theory of gravity [5]. The gravi¬ 
tational waveforms from BH-NS mergers depend on the 
BH mass, the BH spin, the NS mass, and the NS ra¬ 
dius, and thus they will carry information of these bi¬ 
nary parameters |S1[7]. In particular, the information of 
the NS radius could be used for constraining the equa¬ 
tion of state (EOS) and makes a big contribution not 
only to astrophysics but also to nuclear and fundamental 


physics IHHII]. Since the expected signal-to-noise ratio 
in the gravitational-wave detection is not very high, tak¬ 
ing a cross correlation between the observational data 
and theoretical template is an important method to ex¬ 
tract the information from the detected signal at least 
for the next generation gravitational-wave detectors na, 
and hence waveforms for a variety of binary parameters 
have to be derived. 

BH-NS mergers are also proposed as a potential pro¬ 
genitor of short-hard gamma-ray-bursts (sGRBs) [13]. 
If the NS is tidally disrupted during the merger, a hot 
and massive disk with mass > O.OIM© could be formed 
around the remnant spinning BH. This BH-disk system 
could launch a relativistic jet by releasing its gravita¬ 
tional energy through the neutrino emission or an elec¬ 
tromagnetic energy flow in a short time scale ^ 2 s, and 
hence, could be the central engine of the sGRBs. This 
scenario for sGRBs, so-called merger scenario, is suited 
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for explaining its duration and the estimated event rates 
(see Refs. m US] for reviews and references therein). 
For extracting physical information of the central engine 
of sGRB from electromagnetic observations, we have to 
clarify the formation process of the BH-disk system and 
the relation between this system and sGRB. 

A fraction of the NS material would be ejected dur¬ 
ing tidal disruption [H Hlj. Since the NS consists of 
highly neutron-rich matter, r-process nucleosynthesis is 
expected to take place in the ejecta [18], and the emis¬ 
sion powered by decay of the radioactive nuclei would oc¬ 
cur (kilonova/macronova) [TOH^ . This electromagnetic 
counterpart of the binary merger is useful for determin¬ 
ing the position of the source when gravitational waves 
are detected. Also, their light curves will reflect the bi¬ 
nary parameters, and could be useful for extracting the 
physical information of the binary. Furthermore, the r- 
process nuclei produced in the ejecta are considered to 
contribute to the chemical evolution history of the uni¬ 
verse HHIEIES]. 

The property of disks and ejecta, particularly their 
mass, is essential for predicting the electromagnetic coun¬ 
terparts, such as sGRB and kilonova. The mass of the 
matter that remains outside the BH after the merger de¬ 
pends strongly on whether and at which orbital separa¬ 
tion the NS is disrupted by the tidal force of the com¬ 
panion BH [24]. If the tidal disruption occurs at a suf¬ 
ficiently distant orbit from the innermost stable circular 
orbit (ISGO) of the BH, an appreciable amount of mass 
would remain outside the remnant BH. On the other 
hand, if the tidal disruption does not occur or occurs 
near the ISGO of the BH, entire NS material would be 
swallowed by the BH. The orbital separation at which 
the tidal disruption occurs depends on binary parame¬ 
ters, particularly on the mass ratio of the binary, the NS 
compactness, and the BH spin. Thus, we have to clar¬ 
ify the dependence of tidal disruption processes on these 
parameters in order to predict the mass of the remnant 
disk and the ejecta. 

The orbital angular momentum and BH spin can be 
misaligned. A population synthesis study suggests that 
about a half of BH-NS binaries can have an initial con¬ 
figuration in which the misalignment angle between the 
orbital angular momentum and the BH spin is larger than 
45° [25]. A post-Newtonian (PN) study shows that the 
orbit processes due to spin-orbit coupling (or dragging 
of inertial frames) in the presence of misalignment [26] . 
This could affect dynamics of BH-NS binaries not only 
in the inspiral phase but also in the merger phase. 

Recently, a variety of numerical-relativity (NR) sim¬ 
ulations have been performed for BH-NS binaries lilTl 
niiniiiTHiii, and quantitative dependence of the merger 
process on binary parameters have been revealed. Some 
works were done taking into account magnetic fields m 
[371 [38] . nuclear-theory-based EOS [T6[ [TT] [35l [40l [42] . 
and neutrino cooling [40] [42]. However, most of them 
were done for the case that the BH spin is aligned with 
the orbital angular momentum and there are only two 


studies for the misaligned case [3giiT|. Moreover, mass 
ejection has not yet been studied for the misaligned-spin 
case based on NR simulations (there is a study with an 
approximate treatment of general relativity [43]). 

The primary purpose of this paper is to clarify the 
quantitative dependence of the disk formation and the 
mass ejection from BH-NS mergers on BH spin mis¬ 
alignment and nuclear-theory-based EOS by perform¬ 
ing NR simulations systematically. We employed four 
nuclear-theory-based EOSs described by piecewise poly¬ 
trope with four pieces m, and varied the misalignment 
angle between the orbital angular momentum and the 
BH spin, while fixing the NS mass, the magnitude of the 
BH spin, and the mass ratio of the binary as I.35M0, 
0.75, and 5, respectively. 

This paper is organized as follows. In Sec. [n| we de¬ 
scribe methods for computing initial conditions, piece- 
wise polytropic EOS, and models of BH-NS binaries we 
employed in this paper. In Sec. m the formulation and 
the methods of numerical simulations are summarized. In 
Sec. [IV[ we present the numerical results for misaligned- 
spin BH-NS mergers. Einally, summaries and discussions 
of this work are presented in Sec.|V| Throughout this pa¬ 
per, we adopt the geometrical units in which G = c = I, 
where G and c are the gravitational constant and the 
speed of light, respectively. Exceptionally, c is some¬ 
times inserted for clarity when we discuss the velocity 
of ejecta. Our convention of notation for physically im¬ 
portant quantities is summarized in Table |T| The dimen¬ 
sionless spin parameter of the BH, total mass of the sys¬ 
tem at infinite separation, mass ratio, and compactness 


of the NS are defined as y = Q = AfBn/AfNS? 

mo = Mbh + TTns, and C = Mns/Rns, respectively. 
J = {J^^ J'^^ J^) is the total angular momentum of the 
system (see Sec. HIB for its definition). The angle be¬ 
tween the BH spin and the orbital angular momentum, 
Aiit, is defined by 


hilt • COS 


/^Sbh-L\ 

V SbmL J ’ 


( 1 ) 


where L is the orbital angular momentum of the binary, 
which is defined using the total angular momentum, J, 
and the BH spin, Sbh? as L := J —Sbh with Sbb = |Sbh | 
and L = |L|. We use hiit,o to describe itiit at the initial 
condition. Latin and Greek indices denote spatial and 
spacetime components, respectively. 


II. INITIAL CONDITION 
A. Formulation and methods 

We prepare quasiequilibrium states of BH-NS binaries 
as initial conditions of numerical simulations in a similar 
manner to our previous works for aligned-spin binaries 
[SlHIlZ]. Numerical computations are performed by us¬ 
ing multidomain spectral method library LORENE [46] . As 
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TABLE L Our convention of notation for physically important quantities. 


Symbol 

Mirr The irreducible mass of the BH 

Mbh The gravitational mass of the BH in isolation 

Sbh The magnitude of the BH spin angular momentum 

Mns The gravitational mass of the NS in isolation 

i^NS The circumferential radius of the NS in isolation 

mo The total mass of the system at the infinite separation 
Mo The Arnowitt-Deser-Misner mass of the system 

J The total angular momentum of the system 

Q The mass ratio Mbh/Mns 

C The compactness parameter of the NS Mns/^^ns 

X The dimensionless spin parameter of the BH aSbh/Mbh 

Ztiit The misalignment angle of the BH spin 

Ztiit,o ^tiit at the initial condition 


gravitational radiation reaction reduces the orbital eccen¬ 
tricity [48], typical BH-NS binaries settle to a quasicir¬ 
cular orbit. For the case that the orbital separation is 
large, the time scale of radiation reaction, tgw, is much 
longer than the orbital period, Porb? because their ratio 
is given by m 


tgw _ ^ ^ (1 + Qf ( r 
Porb ' 4 Q \6mo J 


( 2 ) 


For the case that the BH has a spin angular momentum 
inclined to the orbital angular momentum, the orbital 
plane precesses due to the spin-orbit coupling effect of 
general relativity, and a closed orbit is not obtained even 
in the absence of radiation reaction [26|. This implies 
that even in one orbital cycle, the gravitational interac¬ 
tion between two objects varies depending on the angle 
between the BH spin and the line connecting two centers 
of mass, and hence, the definition of the quasicircular or¬ 
bit is not trivial. At a large orbital separation, however, 
the precession time scale, Pprec? is also longer than Porb 
as 


-^prec ^ 1 y 7(1 + Q)^ T 
Porb * 4(4(3 + 3) 6mo ’ 

where we take the small-spin limit. Thus, we can neglect 
the orbital precession as well as the gravitational radi¬ 
ation reaction for computing initial data at a large or¬ 
bital separation, and the binary can be regarded approx¬ 
imately as an equilibrium configuration in the comoving 
frame. 

We then compute quasiequilibrium states assuming the 
presence of an instantaneous helical Killing vector field 
with the orbital angular velocity 11, 

e^:=(5,r+ii(5^r- (4) 

For nonspinning or aligned-spin BHs, this reduces to a 
genuine helical Killing vector as far as the radiation re¬ 
action is neglected. Accordingly, the orbital plane can be 


taken to be a plane perpendicular to the rotational axis. 
This does not hold for misaligned-spin BHs, and we do 
not restrict the orbital plane to be a plane perpendicular 
to the rotational axis. Instead, we compute initial data 
by requiring that neither the BH nor NS has the velocity 
component along the rotational axis [36]. This implies 
that the binary is located at the extrema of coordinate 
separation along the rotational axis, where the helical 
symmetry should hold instantaneously. 

The line element in the 3 + 1 form is written as 


ds^ = g^ydx^dx^ 

= —a^dt^ + ^ij (dx'^ + (dx^ + 13^dt) , (5) 


where (a, /3*, and are the lapse function, shift vector, 
and three-diemensional spatial metric, respectively. Ini¬ 
tial data of the gravitational field consist of jij and the 
extrinsic curvature defined by 

Kij •= ( 6 ) 

where is the future-oriented timelike unit normal vec¬ 
tor field to the initial hypersurface. 

We employ the extended conformal thin-sandwich for¬ 
malism |49j |50] with the conformal transverse-traceless 
decomposition of Einstein’s equation m to compute a, 
jij, and Kij. We assume the conformal flatness of 
the spatial metric, jij = the stationar- 

ity of the conformal metric, dt^ij = 0, and the maximal 
slicing condition, K = 0 = dfK. Here, fij and K are 
the flat spatial metric and the trace part of the extrinsic 
curvature, K := Y^Kij^ respectively. To handle a coordi¬ 
nate singularity associated with the BH, we decompose 
the conformal factor 'ijj and a weighted lapse function 
4> := into singular and regular parts using constants 
Mp and as 1^117] 


tA = 1 + 

2rBH 

$ = 1-h rj, 

rBH 


( 7 ) 

(8) 
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where tbr •= |x — xp| is the coordinate distance from 
the puncture located at xp. We also decompose Aij := 
2 p~‘^Kij into regular and singular parts as 

Aij = DiWj + DjWi - pijDkW’^ + Kl, (9) 


where Di denotes a covariant derivative associated with 
fij and the index of is raised/lowered by fij. The 
singular part is given by [52] 


K^. •= _ 

• 2r.2 


XiPf^ + XjPl 


BH 


{fij -XiXj)x^P^ 


BH 


T 3 \fikl^p^ T ^/c] 

^BH 


( 10 ) 


where P := {x^ — Xp)/rBH- The index of P is also 
raised/lowered by fij. The parameters PP^ and Sp are 
constants associated with the linear momentum and spin 
angular momentum of the puncture, respectively (see be¬ 
low). 

The equations to determine 0, and Wi are de¬ 

rived by combining the Hamiltonian constatint, momen¬ 
tum constraint, dtK = 0, and dt^ij = 0 as 

A(f> = -27r'*/)®/9H - (11) 

A/?* + = 167r$V’^i* + 2A^Wj , (12) 

o 

At] = 27r$'*/)^ {pii + 2S) + '^^'tp~^AijA^^, (13) 

O 

AWi + ^bibjW^ = STT-iP^ji. (14) 


where A := p^DiDj. The matter source terms are de¬ 
fined by 

PH :=T^un^n\ (15) 

/ := (16) 

5 *^':= ( 17 ) 


with S := Sij. The asymptotic flatness gives the outer 

boundary conditions as 


<^|oo = ^ioo = r/loo = mu =0. (18) 


other mass parameter, M<^, is determined by the con¬ 
dition that the Arnowitt-Deser-Misner mass and the 
Komar mass agree with each other for stationary and 
asymptotically flat spacetime [58l[59], that is. 



di^dS^ 



difjds^ 


‘2'kMq . 


(19) 


The linear momentum of the puncture, TpH? deter¬ 
mined by the condition that the Arnowitt-Deser-Misner 
linear momentum of the system vanishes as 

P^^ = - J ji^^Sx. (20) 


The spin parameter of the puncture, Sp^ is given in 
Cartesian coordinates: 


Sp — 5'p(sini[ii^, 0, cosi[ji^). (21) 

The magnitude Sp is adjusted for obtaining a desired 
value of BH spin, S'bh, measured on the horizonQ The 
inclination angle is chosen to be 30°, 60°, and 90° in 
this study. Because the direction of orbital angular mo¬ 
mentum does not agree with the axis of helical symmetry, 
alway agree with Aiit ? which is defined as 
the angle between the orbital and BH spin angular mo¬ 
menta. The typical difference is 3°, and we do not adjust 
values of to control itiit in this study. 

The NS matter is assumed to be composed of an ideal 
fluid. The energy-momentum tensor is given by 

:= phUf^u,, + Pg^u, (22) 

where p is the rest-mass density, P is the pressure, 
h = l-\-ePP/p is the specific enthalpy with 5 the spe¬ 
cific internal energy, and is the four-velocity of the 
fluid. The velocity field of the fluid is expected to be ir- 
rotational, because the viscosity of the NS matter is low 
|60l|6T] and the rotational periods of observed NS in com¬ 
pact binaries are not very short (see, e.g., [62|). The zero 
relativistic vorticity condition, or irrotationaliry condi¬ 
tion, is written as 

■= + U^ub (VaUp - VpUa) 

= h~^ - '^iy{hub] 

= 0, (23) 


In contrast to initial data of BH-NS binaries computed 
in the excision framework [531157] , we do not have to give 
inner boundary conditions at the horizon. We also do not 
have to give nonzero boost contributions to the shift vec¬ 
tor [36], because the Arnowitt-Deser-Misner linear mo¬ 
mentum of the system can be set to zero by choosing 
PP^ appropriately. Hence, BH-NS binaries with mis¬ 
aligned spins do not exhibit the center-of-mass motion in 
the puncture framework. 

Free parameters associated with the puncture are de¬ 
termined as follows. The so-called puncture mass, Mp, 
is adjusted for obtaining a desired BH mass, Mbh- The 


where the energy-momentum conservation and adiabac- 
ity condition are used for deriving the second-line expres¬ 
sion (63] [64] . This implies the presence of a velocity po¬ 
tential T such that hu^ = and the elliptic equation 

to determine T is derived from the continuity equation 
V^(pi4^) = 0 together with the helical symmetry. The ir- 
rotational conditions and helically symmetric conditions 


^ Although we distinguish Sp from S'bh, the difference between 
two is at most 0(10“'^). 
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of the specific momentum, C^{hu^) = 0, are combined to 
give 

= —C{= const.), (24) 

which is used to determine h. The integration constant 
C is determined by the condition that the baryon rest 
mass of the NS takes a desired value. As we explain in 
the next section, the specific enthalpy determines all the 
other thermodynamical quantities in the computation of 
initial data. 

The relative location of each component of the binary 
is determined as follows. We fix the binary separation in 
the direction perpendicular to the rotational axis, which 
is chosen to be the 2 ;-axis. The centers of BH and NS 
are put on the xz-plane, but now they are not limited to 
the x-axis0 The orbital angular velocity of the binary 
is determined by the force-balance condition that the NS 
does not move perpendicular to the rotational axis, and 
this amounts to requiring dh/dx = 0 in our coordinates. 
The location of the rotational axis with respect to the bi¬ 
nary components is determined by the condition that the 
magnitude of the orbital angular momentum agrees with 
the value derived by the third-and-a-half post-Newtonian 
formulas for a given value of the orbital angular velocity 
(see Appendix D of [65]). Finally, the binary separation 
along the rotational axis is determined by the condition 
that the NS has no velocity component along the rota¬ 
tional axis, and this amounts to requiring dh/dz = 0 in 
our coordinates [36] . 

B. Piecewise polytropic equations of state 

Since the cooling time scale of NS is much shorter than 
the lifetime of typical compact binaries, we can employ 
a zero-temperature EOS for the NS just before binary 
mergers [66]. Employing a zero-temperature EOS, the 
thermodynamical quantities, such as P, e, and h, can be 
described as functions of p as 

P = P{p) ,£ = £{p) ,h = h{p). (25) 

Erom the first-law of thermodynamics, these quantities 
satisfy relations, 

P 

ds = -izdp, (26) 

p^ 

dh = -dP, (27) 

P 

which determine 5 and h from given P (p), respectively. 


^ We can also fix the locations of both components by forcing them 
to be on the a^-axis. This method respects the instantaneous 
helical symmetry as well as the method adopted in this study. 
Taking the fact that no criteria are available to determine which 
condition gives superior initial data, we simply follow the method 
adopted in [36]. 


In this work, we employ a piecewise polytropic 
EOS [44] to describe a zero-temperature EOS of the 
NS. This is a phenomenologically parametrized EOS, 
which reproduces a zero-temperature nuclear-theory- 
based EOS at high density only with a small number 
of polytropic constants and indices as 

P (p) = Kip^' forpi-i < p < Pi (1 < i < n), (28) 

where n is the number of the pieces used to parametrize 
an EOS. pi is the rest-mass density at a boundary of two 
neighboring i-th and {i + l)-th pieces, Ki is the i-th poly¬ 
tropic constant, and T^ is the Tth adiabatic index. Note, 
here, po = 0 Pn ^ Requiring the continuity 
of the pressure, (1 < — 1), the 

EOS is determined completely by giving pi and T^. 
5 is determined by integrating Eq. ( [^ with the integra¬ 
tion constant, e (0) = 0. It was shown that the piecewise 
polytropic EOS with four pieces reproduces the nuclear- 
theory-based EOSs within ^ 5% errors in pressure for 
the nuclear density range [44] . 

Table [XT] lists the EOSs which we employ in our study. 
We employ the models of the NS EOS which can realize 
the NS with M^s ~ 2Mq which satisfies the recent ob¬ 
servational constraint [671|68]. Eor these models, the NS 
radius is in the range ^11 — 15 km for Mns = I.SSMq, 
which is largely consistent with the recent theoretical and 
observational suggestion [691170]. Eollowing mmm , 
we always fix the parameters of EOS in the subnuclear- 
density region as 

Ti = 1.35692395, (29) 

Ki/c^ = 3.99873692 x 10“® , (30) 

and we set p 2 = g/cm^ and ps = 10^^ g/cm^. Here, 

we insert c for clarity. Instead of giving pi , we give P 2 in 
Table [H] for each EOS, which is the pressure at p = p 2 . 


C. Models 

As we already mentioned, we choose Aiit,o ~ 30°, 60°, 
and 90°. We employ four different piecewise polytropic 
EOSs, APR4, ALE2, H4, and MSI, for each value of 
diit,o- On the other hand, we set the NS mass, Mns, fhe 
mass ratio, Q, and dimensionless spin paramter, y, to be 
fixed values (Mns,Q,x) = (1.35,5,0.75), for which the 
misalignment of the BH spin has a significant effect on 
tidal disruption. Eor all the models, the initial angular 
velocity Oq normalized by the total mass is set to be 
777 , 0 = 0.036. We summarize several key quantities 
for the initial condition in Table mi The label for the 
model denotes the EOS name and the value of itiit,o- 
Specifically, “i30”, “i60”, and “i90” denote the models 
with itiit,o ~ 30°, 60°, and 90°, respectively. Eor all the 
models, we rotate the initial data before we start the 
simulation so that the initial direction of the total angular 
momentum agrees with the direction of the 2 )-axis. 
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TABLE IL The key quantities for piecewise polytropic EOSs M which we employ in this paper. P 2 is the pressure at p = p 2 
shown in the unit of dyne/cm^, T^ is the adiabatic index for each piecewise polytrope, and Mmax is the maximum mass of the 
spherical NS for a given EOS. R 1 . 35 , pi. 35 , and C 1.35 are the radius, the central rest-mass density, the baryon rest mass, 

and the compactness parameter for the NS with Mns = l.SSM©, respectively. 


Model 

logio ^’2 

Ta Tg r4 

Mmax[M0] 

Ri.35[km] 

pi.35[g/cm^] 

M*,i.35[Mo] Cl.35 

APR4 

34.269 

2.830 3.445 3.348 

2.20 

11.1 

8.9x10'^ 

1.50 

0.180 

ALE2 

34.616 

4.070 2.411 1.890 

1.99 

12.4 

6.4x10^^ 

1.49 

0.161 

H4 

34.669 

2.909 2.246 2.144 

2.03 

13.6 

5.5x10^^ 

1.47 

0.147 

MSI 

34.858 

3.224 3.033 1.325 

2.77 

14.4 

4.2x10^^ 

1.46 

0.138 


TABLE III. Key parameters and quantities for the initial con¬ 
ditions adopted in our numerical simulation. The adopted 
EOS, the initial angle between orbital angular momentum and 
the BH spin (ztiit,o), the ADM mass (Mo), and the total angu¬ 
lar momentum (Jo) , respectively. Note that Mns = l.SSM© 
and mo = S.IOM© 


Model 

EOS 

^tiit,on 

Mo[M@] Jo[GMq/c] 

APR4i30 

APR4 

33 

8.04 

63 

APR4i60 

APR4 

63 

8.05 

57 

APR4i90 

APR4 

94 

8.05 

47 

ALF2i30 

ALE2 

33 

8.04 

63 

ALF2i60 

ALE2 

63 

8.05 

57 

ALF2i90 

ALE2 

94 

8.05 

47 

H4i30 

H4 

33 

8.04 

63 

H4i60 

H4 

63 

8.05 

57 

H4i90 

H4 

94 

8.05 

47 

MSli30 

MSI 

32 

8.04 

63 

MSli60 

MSI 

63 

8.05 

57 

MSli90 

MSI 

93 

8.05 

48 


tion equations of the constraints have a non-propagating 
mode. In the Z4c formulation, by introducing a new vari¬ 
able, 0 , the evolution equation of the constraints changes 
entirely to a wave equation. Therefore, the violation of 
the constraints propagates away, and the local accumu¬ 
lation of the numerical error can be suppressed. 

In the Z4c formulation, we evolve the conformal factor, 
W := 7 “^/^, the conformal three-metric, jij = 
a variable slightly modified from the trace of the extrinsic 
curvature, K := JC — 20, the conformal trace-free part 
of the extrinsic curvature, Aij = 7 “^/^ [Kij — 
an auxiliary variable, r% and the new variable, 0. Here, 
7 := The evolution equations are written as. 


a(k + 20 ) - 

{dt - 13'"dk) % = -2aAij + jikdj/S'" + jjkdi/S'" 


{dt-l3'di)W=-W 


^lijdkl3'", 


(31) 

(32) 


III. METHODS OF SIMULATIONS 

Numerical simulations are performed using an adap¬ 
tive mesh refinement (AMR) code SACRA [72]. Here, we 
employ a Baumgarte-Shapiro-Shibata-Nakamura (BSSN) 
formulation partially incorporating Z4c prescription m- 
The gauge conditions, the numerical scheme, and the 
diagnostics are essentially the same as those described 
in inilT]. 


{dt - /3'di) K = -D'Dia 


i(A- + 26)^ 


+ 47r (pH + S)] 

+ aKi (1 - K 2 ) 0, 


{dt - 13^dj) 0 = 


it, {if- 


AuA"^ + - 

-r 3 


- a {Stt/Oh + (2 + ^ 2 ) 0} 


(33) 

(i^ + 20)^1 

(--/r-o(^g4) 


A. Formulation and Numerical Methods 

We solve Einstein’s evolution equation partially incor¬ 
porating the Z4c formulation (see m for our prescrip¬ 
tion) with a moving puncture gauge. The Z4c formula¬ 
tion is a modified version of the BSSN-puncture formu¬ 
lation, introducing a new variable. In the BSSN formu¬ 
lation, the long-term simulation causes a slow accumula¬ 
tion of the numerical error that leads to gradual violation 
of the constraints, which should be zero if the evolution 
equations are solved exactly. This accumulation comes 
from the fact that in the BSSN formulation, the evolu- 


{dt - 13'"dk) Aij = -W^ 


D^DjOL Dk^ 


+ 

+ 


W a [ Rij - -^ijR 


a 


(x + 20) 


Aii - 2AikA'" 


+ ^ikdjp^ + Ajkdip^ 

- lAijdkf3'", 


f Sij - - 


( 35 ) 
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{dt - p^dj) r = -2A^^dja + 2a 

- ( 2 k + 0 ) - 87rf 


+ + If^djdkP'^ - f^afe/3* 

+ - 2aKi (P - f^) , (36) 


where Di denotes a covariant derivative associated with 
rjj = —and /^i and k ,2 are coefficients asso¬ 
ciated with the constraint damping. An overall factor, 
g-(r/ro) ^ multiplied in the right-hand side of Eq. (34) 
so that 0 plays a role only in the inner region of the 
simulation box. In our simulation, we set = /t :2 = 0, 
and ro = A/2, where L is the size of the computational 
domain on one side (see Table. IV). The spatial deriva¬ 
tives in the evolution equations are evaluated by fourth- 
order centered finite differencing except for the advection 
terms, which are evaluated by fourth-order upwind finite 
differencing. A fourth-order Runge-Kutta scheme is em¬ 
ployed for the time evolution. 

Following we employ a moving-puncture gauge in 
the form 


{dt - P^di) a = -2aK, (37) 

{dt - p^dj) k = \b\ (38) 

{dt - p^dj) = {dt - p^dj) r - nBB\ (39) 


where is an auxiliary variable, and 773 is a coefficient 
introduced to suppress a strong oscillation of the shift 
vector. In this work, we set 77 B = O.IG/Mq. 

The EOS is basically the same as those described in [6l 
[ 7 ]: We decompose the pressure and the specific internal 
energy into a zero-temperature part and a thermal part 
as 


P — Pcold Pth 5 ^ — ^cold “ 1 “ '^th* 


atmosphere density is set to be patm = 10 ^^Pmax ^ 
10^ g/cm^ for the inner region, r < Rent, and patm = 
10“^^Pmax (Acrit/^)^ fol* 1^^ OUtcr rCgioU, r > Rcrit, with 
Rcrit ~ A/16. The total rest mass of the atmosphere 
is always less than 10 “^Mq, and hence we can safely 
neglect the effect of the artificial atmosphere as far as 
appreciable tidal disruption occurs. 


B. Diagnostics 


We estimate the mass and the spin angular momentum 
of BH with misaligned spin assuming that the deviation 
from Kerr spacetime is negligible in the vicinity of the 
BH, at least for the case that the separation of the binary 
is large or the system is approximately regarded as a 
steady state. In a steady state with stationary slicing, the 
event horizon (EH) agrees with apparent horizon (AH). 
Thus, the irreducible mass of the BH is determined by 
the area of the AH, Aah, as 


Mirr 



In Kerr spacetime, a relation 

2Mi (1 - yr^) 


MgH — 




(42) 


(43) 


holds between the gravitational and irreducible masses 
of the BH. If the irreducible mass and the BH spin are 
known, we can calculate the BH mass using this relation. 
To determine the spin angular momentum, we use the re¬ 
lation between the spin and the intrinsic scalar curvature 
of the horizon, following [78]. In Kerr spacetime, 

is written as 


(2) ... ^ 2 (7 + x") (7 - Sx^cos^^) 

-^BH (7 + X^COs26>)^ 


(44) 


Here, Acoid and Scoid are functions determined by the 
piecewise polytropic EOS. Then, the thermal part of the 
specific internal energy is calculated by eth = s — ^coid, 
where 5 is given from the hydrodynamics. Finally, we 
determine the thermal part of the pressure using a simple 
T-law, ideal-gas EOS as 

Ath = (Eth - 1) (41) 

where Tth is an adiabatic index for the thermal part, 
which we choose Tth = 1-8. As is discussed in Appendix 
A. of [17], the difference of the disk mass and the ejecta 
mass among the different values of Tth is expected to 
be small compared to the numerical errors due to finite 
gridding, thus, we do not study the dependence on Tth 
in this paper. 

Since the vacuum is not allowed in any conservative 
hydrodynamics scheme (see [72] for details), we put an 
artificial atmosphere of small density outside the NS. The 


Here, r+ = 1 + a/ 1 — is a normalized radius of the 
EH, and 0 is the latitude in Boyer-Lindquist coordinates. 
The minimum and maximum values of A^j^ at ^ = 0 
and 7r/2, i.e., at the pole and the equator of the BH are, 
respectively. 


d( 2 ) 

^min 




-1 + 2 ^ 1 -x^ 

2M? 

irr 

2 (-2 + x' + 2 ^ 1 -x^) 


(45) 

(46) 


Solving these equations with respect to x, we have 


Amin 


— 1 “ 1 2 ’ 


X^ = 

A max 


-2 + 2^2M?^R^l 


(47) 

(48) 


















TABLE IV. Setups of the grid structure for the simulation 
with our AMR algorithm. Ax is the grid spacing at the finest- 
resolution domain. Rdiam is the semimajor diameter of the NS 
in the direction perpendicular to the axis of helical symmetry. 
L is a half of the edge length of the largest domain. Ao = tt/Qq 
is the gravitational wavelength of the initial configuration. 


Model 

Ax/Mo 

Rdiam/Ax L j Xq L (km) 

APR4i30 

0.0134 

101 

2.357 

2444 

APR4i60 

0.0133 

102 

2.343 

2429 

APR4i90 

0.0136 

100 

2.398 

2486 

ALF2i30 

0.0153 

102 

2.687 

2786 

ALF2i60 

0.0153 

102 

2.687 

2786 

ALF2i90 

0.0156 

101 

2.729 

2829 

H4i30 

0.0172 

102 

3.032 

3144 

H4i60 

0.0172 

102 

3.032 

3144 

H4i90 

0.0173 

101 

3.046 

3158 

MSli30 

0.0186 

102 

3.273 

3394 

MSli60 

0.0188 

101 

3.308 

3430 

MSli90 

0.0188 

101 

3.308 

3430 


Using these relations, we can estimate the BH spin ap¬ 
proximately from the intrinsic scalar curvature of the AH. 
The direction of the BH spin can be determined from the 
location of We note that the direction of the BH 

spin is not gauge invariant in this definition. However, 
we expect that it gives a reasonable measure of the spin 
direction if the tidal forces are negligible [TEllTg]. For all 
the models, the values of Xmin and Xmax agree with each 
other up to significant 3) digits . Thus, we use the 
value of Xmax for calculating Mbh and Sbh in this paper. 

We define the total angular momentum of the system, 
J = (J^, J^, J^), from a rotational invariance of the grav¬ 
itational Hamiltonian at spatial infinity as [45], 

X := {Kl - Kji) dSi, (49) 

OTT 

where is the spatial Levi-Civita tensor. 

Finally, we define the time at which the binary merges. 
For this purpose, we define the rest mass inside the AH 
as 

^<AH •= / (50) 

Jr<rAH 

where, tah = ^ah (^5 is the radius of the AH as a func¬ 
tion of the angular coordinates, and p* = pau^^^ is the 
conserved rest-mass density. Then, we define the merger 
time, t merger ^s thc time at which M<ah ^ 10 “^ Mq is 
achieved. 


C. Setups for AMR grids 

In SACRA, the Einstein and hydro dynamical equations 
are solved in an AMR algorithm described in [72|. Here, 
we briefly describe the settings for AMR grids, and the 
details are found in m- In this work, we prepare 


nine refinement levels with different grid resolutions and 
domain sizes. Each domain is composed of the uni¬ 
form vertex-centered cubic grid with the grid number 
{2N -h 1, 2N 1, 2N + 1) for (x, p, z). We always chose 
N = 60 for the best resolved runs in this work. We 
also performed simulations with N = 40 and 48 to check 
the convergence of the result. As described in [72], the 
AMR domains are classified into two categories: One is 
the coarser domains which cover wider regions with their 
origin fixed approximately at the center of the mass of the 
system. The other is the finer domains which cover the 
regions around the BH or the NS and comove with it. We 
set four coarser domains and five pairs of finer domains 
for all the simulations which we performed in this paper. 
The grid spacing for each domain is hi = I//(2W), where 
2L is the edge length of the largest cubic domain and I 
is the depth of the domain. 

Table [rV] summarizes the parameters of the grid struc¬ 
ture for the simulations. In all the simulations, the semi¬ 
major diameter of the NS in the direction perpendicular 
to the axis of helical symmetry is covered with « 100 grid 
points. For A = 60, the total memory required for the 
simulation with 14 domains is about 35 GBytes. We per¬ 
form all the simulations using personal computers of 64- 
128 GBytes memory and 6-24 processors with OpenMP 
library. The typical computational time required to per¬ 
form one simulation is 9 weeks for the 24 processors case. 


IV. NUMERICAL RESULTS 

We present numerical results of our simulations in this 
section. 


A. Orbital Evolution 

Figure plots bird’s-eye views for the evolution of 
the coordinate separation Xg^p := xj^g — Xgj^ for mod¬ 
els H4i30, H4i60, and H4i90. Here xj^g denotes the lo¬ 
cation of the maximum rest-mass density, and Xgj^ is 
the location of the puncture. This figure shows that the 
number of the orbits before the merger decreases with 
the increase of itiit,o as 8.5, 7.5, and 5.5, respectively 
for models H4i30, H4i60, and H4i90. This dependence 
stems primarily from the general relativistic spin-orbit 
interaction, which is well known for inducing an “orbital 
hung up” effect [261I75H77] . The additional energy of the 
spin-orbit interaction is written, in the leading order, as 

^so = 4^L.Sbh, (51) 

r'^ Mbh 

and hence the spin-orbit interaction essentially weak¬ 
ens the attractive force of gravity if L • Sbh > 0- For 
this situation, the angular velocity of the binary, U, 
decreases, and so does the luminosity of gravitational 
waves, which is proportional to The reduction of 
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H4i30 


H4i60 


H4i90 





FIG. 1. Evolution of the orbital separation Xgep := ~ ^bh of binaries with (Q,Mns,x) = (5, l.SSM©, 0.75). The left, 

middle, and right panels show the results with H4i30, H4i60, and H4i90, respectively. 
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FIG. 2. Evolution of itiit for the models with H4. 


FIG. 3. Evolution of the precession angle, 0prec, for the mod¬ 
els with H4. 


the gravitational-wave luminosity makes the approaching 
velocity smaller, and thus the time to merger becomes 
longer. Since L • Sbh cx: cositiit, this effect can be signifi¬ 
cant when the BH spin is aligned with the orbital angular 
momentum, and thus, the binary with a small value of 
hiit,o merges later. 

The figure also illustrates that the orbits of the bina¬ 
ries are precessing. This is also primarily due to the spin- 
orbit interaction. For models H4i30, H4i60, and H4i90, 
the elevation angles of the orbits measured from the xy- 
plane are always ~ 15°, 30°, and 45°, respectively. We 
note that gravitational waves are radiated primarily to 
the direction of the orbital angular momentum, and the 
direction of the total angular momentum J changes dur¬ 
ing the inspiral phase due to the gravitational radiation 
reaction. However, the angle between J and the 2 ;-axis is 
always smaller than 5°. 

It is known that, at least approximately, itiit is a con¬ 
stant of motion [26] for the case that the spin of the NS 
is absent. This feature is also seen in our simulation: 
see Fig. in which we plot the time evolution of itiit- 


Irrespective of itiit,o and EOS, we indeed find that itiit 
approximately keeps their initial values and their fluctu¬ 
ation is smaller than ^ 3° irrespective of models. There¬ 
fore itiit would be regarded approximately as the value 
which is determined when the BH-NS binary was born, 
even just before the merger. Also this property ensures 
that setting of the simulation models is well-defined; the 
models with different values of Aiit,o describes entirely 
different physical systems. 

Next, we analyze the evolution of the precession angle. 
We define the precession angle of the orbit as 


0prec •— COS 


(L X Sbh) • y 
|L X Sbh| 


(52) 


Here, y is the coordinate basis of the ^-axis. The evolu¬ 
tion of 0prec is shown in Fig. up to the time of merger. 
We also plot the value obtained by integrating the leading 
post-Newtonian formula of the precession angular veloc- 
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ity US], 

<ec = ^(3Q-'+4), (53) 

using instantaneous values of J and r for the simulation 
of the model H4i60. 

For models H4i30, H4i60 and H4i90, the final values 
of 0prec are ~ 290°, 235°, and 145°, respectively, while 
the number of the orbits are ^ 8.5-5.5. We find that the 
precession angular velocity cjprec computed by the time 
derivative of 0prec is always smaller than orbital angu¬ 
lar velocity by an order of magnitude. The evolution of 
0prec agrees quantitatively with the one calculated with 
the leading-order post-Newtonian formula, despite the 
difference in the gauge condition. Figure shows that 
the final value of 0prec is smaller for the case that diit,o 
is larger. It is simply because a longer inspiral phase is 
realized for binaries of a smaller value of itiit,o- 

B. Tidal disruption 


pointed out in mi that this may prevent the collision of 
the tidal tail. However, the tidal tail still collides with 
itself to form a disk or torus for this model. This might 
be due to the difference of the initial parameters of the 
binary, that a binary with larger mass ratio and dimen¬ 
sionless spin parameter are employed in mi, for which 
larger elevation angle is expected to be achieved. More¬ 
over, we performed the simulation for a longer time after 
the merger than in m, and this also makes more chances 
for the tidal tail to collide with itself. The disk appears 
to be misaligned with the BH spin by ~ 20° at least at 
« 10 ms after the tidal disruption, while the BH spin is 
aligned approximately with the 2 )-axis. 

For the model H4i90, the NS is tidally disrupted very 
weakly by the BH and a tiny tidal tail is formed. Since 
the tidal disruption occurs at a close orbit to the ISCO 
or perhaps inside the ISCO of the BH, most of the NS 
material is swallowed by the BH. Thus, only a tiny ac¬ 
cretion disk is formed and the total amount of the ejecta 
is not appreciable for this model. 

Figure plots the time evolution of the total rest mass 
outside the BH defined by 


Figure shows the rest-mass density together with 
the location of apparent horizons for models H4i30 and 
H4i60. The images are generated using a volume ren¬ 
dering method truncating the density below 10^^ g/cm^, 
and the color on the AH surfaces descri bes th e value of 
2-dimensional Ricci scalar on it (cf. Sec. HIB). 

For the model H4i30, the NS is tidally disrupted form¬ 
ing a one-armed tidal tail around the BH. An efficient 
angular momentum transport process works, and a frac¬ 
tion of the NS material becomes gravitationally unbound 
during this phase. At a few milliseconds after the onset 
of the tidal disruption, a large fraction of the NS material 
is swallowed by the BH, and ~ 20% of the NS material 
remains outside the BH. The inner part of the tidal tail is 
subsequently wound around the BH, and a disk with its 
radius ~ 150 km is formed. Also, some material, which 
was not able to get enough kinetic energy to escape from 
the system, falls back to the disk continuously. 

Initially, the tidal tail around the BH is tilted with 
the angle less than 15°, which reflects the elevation 
angle of the orbit. Also, the tidal tail is slightly warp¬ 
ing due to the spin-orbit interaction. However, the disk 
formed finally is nearly aligned with the BH spin, and the 
morphology of the disk for this model resembles the disk 
formed with aligned-spin BH-NS mergers. The remnant 
BH spin axis agrees approximately with the direction of 
the initial total angular momentum, i.e., the ^-axis. 

For the model H4i60, general features of the tidal dis¬ 
ruption of the NS is similar to the model H4i30, while 
the precession of the tidal tail is more appreciable for 
this model. More than 90% of the NS material falls 
into the BH in a few milliseconds after the tidal disrup¬ 
tion for this model. This is a result of the fact that 
the tidal disruption occurs in the vicinity of the ISCO of 
the BH. The elevation angle of the tidal tail measured 
from the x^-plane is different for each part, and it is 


Af>AH := / p*d^x, (54) 

Jr>rAii 

for a variety of EOSs and itiit,o- From the comparison 
among the models with the same value of Aiit ,0 5 we find 
that the values of M>ah increases in the order of APR4, 
ALF2, H4, and MSI, i.e., in the order of the compactness. 
This result shows that the tidal disruption occurs at a 
more distant orbit for the case that the compactness of 
the NS is small. This dependence of the tidal disruption 
on the compactness is the same as the dependence which 
was found in the study on BH-NS mergers with aligned 
spins [6l[7]. 

For a fixed EOS, the values of M>ah after the merger 
decrease as the values of Aiit,o increase, and thus the 
increase of Aiit,o prevents the tidal disruption of the 
NS. This result agrees qualitatively with the result of 
previous study on the misaligned-spin BH-NS binary 
merger Esum, and is explained primarily by the reduc¬ 
tion of the spin-orbit coupling. According to the study on 
the aligned-spin BH-NS merger [7], a larger mass remains 
outside the BH after the merger for the case that the BH 
spin is larger and parallel with the orbital angular mo¬ 
mentum. This is because the spin-orbit interaction works 
as a repulsive force and the ISCO radius of the BH be¬ 
comes small for this situation. Since the spin-orbit inter¬ 
action energy is proportional to L • Sbh = AS'ERCOshut^o 
at the leading PN order, the spin-orbit interaction is 
weakened as itiit,o increases. Thus, the increase of Aiit,o 
enlarges the ISCO radius of the BH effectively and re¬ 
duces the remnant mass after the merger. 

We summarize the value of M>ah at ~ 10 ms after 
the onset of merger in Table [V| We compare these results 
with the fitting formula for the aligned-spin case obtained 
in [80]. Since the spin-orbit interaction is proportional 
primarily to L • Sbh, we compare our numerical results 
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FIG. 4. Snapshots of the volume rendered density map as well as the location of AHs at ~ 5 ms after the onset of merger for 
the models H4i30 (left panel) and H4i60 (right panel). The inner regions of the simulation with side lengths of ~ 300 km are 
shown. 




FIG. 5. Evolution of the rest mass outside the apparent horizon M>ah for the models with Ztiit,o ~ 60° (left figure) and with 
H4 (right figure). The result for the aligned-spin case with Q = 5 are picked up from [T7]. 


with the values derived by the fittirg formula using the 
effective spin parameter defined by|j 

Xeff = xcositiit,o- (55) 

The deviations of the value calculated by fitting formula 
Afit = |Mfit - M>AH|/^fit are within 50% for M>ah ^ 
O.IMq and within 30% for M>ah ^ 0.2Mq. As far as 
the value of M>ah is larger than O.IMq, M^t gives a 
reasonable estimate for M>ah- 


^ In |41l 1811182| , the effective spin parameter is defined by a differ¬ 
ent form. 


C. Disk formation and mass ejection 

1. Mass of ejecta and disk 

Next we evaluate the mass of disk and ejecta, which 
are the key quantities for the electromagnetic emission 
from the remnant of BH-NS mergers. We calculate the 
ejecta mass by 

Meje := / p*(fx. (56) 

J r>rAU,Ut< — l 

Here, we assume that the contribution of the internal en¬ 
ergy of the ejecta is negligible for estimating the unbound 
material. Then, we define the rest mass of remnant disks 
by 

Mdisk ■= Af>AH - Mej. 


(57) 
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TABLE V. The list of M>ah, Mdisk, Mgje, '^ave, and Ve^e- The subscripts 10 ms and 20 ms denote the values evaluated at 10 ms 
and 20 ms after the onset of merger, respectively. The results for the aligned-spin case are obtained from HU. - implies that 
we were not able to take the data for them. For the model APR4iO, the simulation was stopped before t — tmerge ~ 20 ms. For 
the model APR4i90, the mass of the disk and ejecta are so small that accurate values cannot be derived for them. For the i90 
models, the data for i was not output. 


-Tf>AH, lOms [TfQ] Afdisk, 10ms [TIq] Mdisk,20ms[T/'Q] Meje ,10ms [TIq] '?^ave,10ms [c 

-:————-:— _ ___—- n -:———- 


Model 


'^eje,10ms [cj 

0.099 

0.057 

0.078 


8 X 10“^ 

5 X 10“^ 
1 X 10“^ 
< 10 “^ 


APR4iO 

APR4i30 

APR4i60 

APR4i90 


0.068 
0.022 
4 X 10“ 


0.24 
0.16 
0.026 
2 X 10“^ 


0.059 
0.017 
2 X 10“' 


0.014 
2 X 10“ 


0.26 

0.30 

0.27 

0.24 


ALF2iO 

ALF2i30 

ALF2i60 

ALF2i90 


0.20 
0.13 
0.016 
1 X 10“' 


0.16 
0.10 
0.013 
< 10 


-4 


0.046 
0.033 
0.010 
< 10 


-4 


0.21 

0.27 

0.28 

0.26 


0.15 

0.17 

0.048 


H4i0 

H4i30 

H4i60 

H4i90 


MSliO 

MSli30 

MSli60 

MSli90 


0.32 
0.25 
0.084 
3 X 10“' 


0.27 
0.21 
0.072 
2 X 10“' 


0.21 
0.19 
0.061 
2 X 10“ 


0.050 
0.042 
0.012 
1 X 10“ 


0.22 

0.27 

0.25 

0.28 


0.18 

0.21 

0.14 


0.36 

0.30 

0.18 

0.022 


0.28 

0.23 

0.14 

0.012 


0.22 

0.19 

0.11 

0.011 


0.079 

0.070 

0.041 

0.010 


0.24 

0.28 

0.27 

0.27 


0.19 

0.23 

0.21 


^diskt^sun] 


^ejet^sup] 


0 

0 

&_ 
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"O 
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FIG. 6. The contour for Mdisk (left panel) and Meje (right panel) evaluated at « 10 ms after the onset of merger in the plane 
of the NS compactness C and initial value of Ztiit- 


We note that Mdisk is described as Mbd in m- 

We list the values of Mdisk and Meje at 10 ms and 
^ 20 ms after the onset of merger in Table [V| This shows 
that Mdisk and Meje monotonically decrease with the in¬ 
crease of itiit,o- This reflects the fact that the effective 
IS CO radius of the BH increases with the increase of 
^tiit,o- We find that Mdisk and Mgje with ^ 30° are 
appreciably smaller than those for itiit,o = 0°. In partic¬ 


ular, little amount of disk and ejecta are produced for all 
the EOS with itiit,o ~ 90°. For the moderate misalign¬ 
ment angle, itiit,o ~ 60°, the values of Mdisk and Meje 
are sensitive to the EOS: For the MSI EOS, the disk 
with Mdisk > O.IMq and ejecta with Meje > 1 O“^M 0 
are produced. On the other hand, Mdisk < 10 “^Mq and 
Meje < 1 O“^M 0 for the APR4 EOS. 

To clarify the dependence of Mdisk and Meje on itiit,o 
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FIG. 7. Mdisk (left panel) and Mgje (right panel) evaluated at ~ 10 ms after the onset of merger as a function of an effective 
spin parameter Xeff = xcosztiit,o- The plotted lines are the linear interpolation for the results of misaligned BH-NS mergers 
obtained in this paper, and the plotted points are the results for aligned-spin BH-NS mergers with Q = 5 obtained in HU. 
“a5” and “a375” denote the aligned models with x = 0-3 and x = 0.375, respectively. 


and EOS, we plot contours of Mdisk and Meje at ~ 10 ms 
after the onset of merger as functions of itiit,o and the 
compactness parameter C of the NS in Fig. The depen¬ 
dence of Mdisk and Meje is clear: Both of them decrease 
monotonically with the increases of C and itiit,o- For a 
moderate value of compactness C = 0.160, itiit,o should 
be smaller than 50° for Mdisk to be larger than O.IMq. 
On the other hand, Mdisk ^ O.IMq even if itiit,o ~ 70° for 
a stiff EOS that realizes C = 0.140. For a soft EOS with 
which C > 0.175, Mdisk is smaller than O.IMq for any 
value of itiit,o- Meje > O.OIMq is possible for itiit,o < 85° 
with C = 0.140, for itiit,o < 65° with C = 0.160, and for 
^tiit,o < 30° with C = 0.175. 

Figure [^compares Mdisk and Meje obtained by numeri¬ 
cal simulations for aligned-spin BH-NS mergers m with 
those for the misaligned-spin cases. Each line describes 
the results of Mdisk and Meje for the misaligned-spin BH- 
NS mergers interpolated linearly for Xeff- Each point in 
Fig. shows the results of the aligned-spin BH-NS merg¬ 
ers with the same mass ratio {Q = 5) and the same EOS 
as we employed in this paper, but with smaller BH spin 
X = 0.5. We also plot a new result for the model with 
(5 = 5, H4 EOS, and x = 0.375. For both Mdisk and 
Meje, the results of the aligned-spin case agree approx¬ 
imately with the interpolated line in the error margin 
due to the finite grid resolution (see Appendix , while 
slightly larger mass is realized for the misaligned-spin 
case. The slope between H4a375 and H4a5 also agrees 
with the slope of the interpolated line for H4. Excep¬ 
tionally for models with ALF2, results obtained by the 
aligned-spin BH-NS mergers deviate from the interpo¬ 
lated plot lines by ^ 100% for Mdisk and « 30% for 
Meje, which might be a little bit larger than the error 
margin due to the finite grid resolution. We note that 
the deviation depends on the interpolation method and 
a systematic error associated with employing Xeff- 


2. Disk morphology and accretion 

Figure shows snapshots of the density profiles at 
20 ms after the onset of merger for selected models 
(see figure caption). We find that the size of the re¬ 
gion with p > 10^ g/cm^ is always ^ 150 km, while 
that with p > 10^^ g/cm^ depends on the mass of 
the disk: For the larger disk mass, the size of the re¬ 
gion with p > 10^^ g/cm^ becomes larger. The maxi¬ 
mum density of the disk also becomes larger for larger 
disk mass: For instance, the maximum density exceeds 
10^^ g/cm^ for the model H4i30, while there is no region 
with p > 10^^ g/cm^ for the model H4i60. The ratio of 
the disk height to disk radius depends only weakly on the 
binary parameters, and it is always ^ 0.3-0.5. 

For the models with Aiit,o ~ 30° ^ iFe rotational axes 
of the disks are approximately aligned with the ^-axis. 
However, the density distribution on the xy cross sec¬ 
tion is not axisymmetric and an approximately station¬ 
ary spiral-shape shock wave is seen in the disks. On the 
other hand, for itiit,o ~ 60°, the disks are misaligned with 
the 2 :-axis. Also, in a similar manner to i30 models, the 
disks have a non-axisymmetric structure for itiit,o ~ 60°. 

To quantify the misaligned structure of the disk, we 
define the total angular momentum of the disk as 

■^disk = / pX^'^XjUkd^x, (58) 

J r>rAH,tit> —1 

where Ui = hui^ and we plot the evolution of the tilt angle 
of measured from the direction of the BH spin after 
the merger in Fig[^ This shows that the tilt angle of J^isk 
at ~ 5 ms after the onset of merger is ~ 15° and ~ 30° for 
models with Aiit,o ~ 30° and itiit,o ~ 60°, respectively. 
These tilt angles of reflect the elevation angle of 

the orbits just before the merger and gradually decrease 
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FIG. 8. The density profiles of the accretion disk at 20 ms after the onset of merger for selected models. The left, middle, 
and right columns show the plots for the xy, xz, and yz-planes, respectively. The top, second top, third top, and bottom raws 
show the plots for models ALF2i30, H4i30, H4i60, and MS 1160, respectively. 
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FIG. 9. Evolution of the tilt angle between fhe BH 

spin for the models with Ztiit,o ~ 30° and 60°. 
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FIG. 10. The density profiles on the yz cross section for the 
model MS 1130 at ~ 40ms after the onset of merger. 


as the system evolves. Figure shows that these tilt 
angles appear to be slightly larger than those expected 
from Fig. 1^ We should note that the tilt angle shown in 
Fig. does not really describe representative values for 
the tilt angle of the dense part of the disk but indicates 
the tilt angle of the remnant matter including the tidal 
tail with larger orbital radii (see below). 

Figure 10 shows a larger-scale density profile 
1500 km) on the ^ 2 ;-plane for the model MSli60 at 

40 ms after the onset of merger. While there is a 
relatively dense torus in the central region 200 km), 
the tidal tail is widely spreading with an elevation angle 
« 30°. This tidal tail is not unbound although it has a 
large orbital angular momentum, and thus it contributes 
to the tilt of which we defined above. 


Figure pr] shows the density profiles of the disk on the 
yz cross section for the model MSli60 at ~ 10,40, and 
70 ms after the onset of merger. While the dense part 


with p > 10^ g/cm^ is tilted by ^ 30° from the xy- 
plane at 10 ms after the onset of merger, the tilt angle 
of the disk decreases gradually as the system evolves, 
and its axis is approximately aligned with the 2 ;-axis at 
70 ms after the onset of merger. The same feature 
can also be seen on the plot of the x 2 ;-plane, and thus 
we conclude that the disk has a tendency to align with 
the BH spin during the evolution. The time scale for 
the disk to align is 50 ms and is comparable to or 
even bit shorter than the time scale of the disk preces¬ 
sion tprec ^ 100 Ills (rdisk/150 km)^. In the presence of 
fluid viscosity, the so-called Bardeen-Petterson effect IHZl 
is known as a mechanism for the disk to be aligned with 
the BH spin. Since any effect of viscosity is not taken 
account, the Bardeen-Petterson effect cannot play a role 
in our simulation. However, we suspect that Bardeen- 
Petterson-like effect induced by a purely hydrodynami- 
cal mechanism, such as angular momentum redistribu¬ 
tion due to a shock wave excited in a non-axisymmetric 
manner of the disk, should work in the disk for the align¬ 
ment. To summarize, a disk with tilt angle of ~ 20°-30° 
would be formed for models with itiit,o ~ 60°. However, 
the dense part of the disk is subsequently aligned with 
the BH spin in ^ 50 ms while the tidal tail with large 
orbital radii keeps its elevation angle. 

Figurefl^plots the time evolution of M^isk foi* the mod¬ 
els with the MSI EOS and itiit,o ~ 30°. This figure shows 
that Mdisk gradually decreases after a steep initial de¬ 
crease. This shows that the infall of the material into 
the BH continues for a long time scale. This is induced 
primarily by a hydrodynamical process associated with 
the non-axisymmetric torque in the disk, as we can infer 
from Fig. Also, the fallback of the matter could give 
an impact to the disk material. To check that numerical 
viscosity is not a main source of the angular momentum 
redistribution, we plot the mass accretion rate of the disk 
material into the BH for three different grid resolutions 
in Fig. Here the accretion rate is defined as the time 
derivative of M<ah- We find that the value of M<ah de¬ 
pends only very weakly on the grid resolution. Since the 
numerical viscosity should depend on the grid resolution, 
this result shows that the contribution of the numerical 
viscosity to the mass accretion is negligible, and we can 
safely consider that the accretion is induced by a physical 
process. The non-axisymmetric and dynamical feature of 
the disk would be responsible for this. 

Lower panels of Fig. plot the time scale of the accre¬ 
tion, tacc •= Ardisk/^<AH, for the models with the MSI 
EOS and for qiit,o ~ 30 °. We find that it depends only 
weakly on Aiit,o and EOS and tacc ~ 100ms for all the 
models at 10-20 ms after the onset of merger. The time 
scale increases as the system relaxes, but still, tacc is as 
short as ~ 200 ms at 50 ms after the onset of merger. 
In other words, the mass accretion rate is as large as 
^ 0.5-1Mq/s for Mdisk = O.IMq even at 50 ms after 
the onset of merger. This shows the importance of the 
non-axisymmetric structure of the disk that is preserved 
for a long time scale and governs the angular momentum 
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FIG. 11. The density profiles on the yz cross section for the model MSliGO at 10,40, and 70 ms after the onset of merger. 
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FIG. 12. Evolution of Mdisk and the accretion time scale tacc = Mdisk/M<AH for the models with MSI (left figure) and with 
^tiit,o ~ 30° (right figure). 


transport. 


velocity of the ejecta as, 


3. Ejecta morphology and velocity 



(59) 


In Fig.[^ we plot the image of the ejecta with volume 
rendering for models H4i30, H4i60, and ALF2i60. We 
find that for the case that the ejecta is massive (H4i30), 
it has a crescent-like shape with its opening angle ~ 180°. 
On the other hand, for the case that the ejecta has small 
mass (ALF2i60), the opening angle becomes larger than 
360° winding around the center of mass. We also find 
that the ejecta is warped due to the orbital precession 
for this case. In particular, the rear-end collision for the 
ejecta with its opening angle larger than 360° is less pro¬ 
nounced in the misaligned-spin case than in the aligned- 
spin case m- 

After the matter becomes gravitationally unbound, the 
ejecta expand in an approximately self-similar manner. 
To analyze the feature of the ejecta, we define an average 


where Teje is kinetic energy of the ejecta, defined follow¬ 
ing [83]. We also compute the linear momentum of the 
ejecta defined by 



and calculate its magnitude by 


P ■ — 

^eje — 



(60) 


(61) 


Using Pejei we define the bulk velocity of the ejecta Veje 

by 


_ -Peje 

~ M ■ ' 


( 62 ) 
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FIG. 13. Evolution of the accretion rate, M<ah, to the AH 
for MSliSO but with different grid resolutions. 


Here, note that may reflect the morphology of the 
ejecta: The linear momentum of the ejecta vanishes if 
its morphology is isotropic, while the value of 'r’eje be¬ 
comes close to i;ave if the mass is ejected coherently to a 
particular direction. 

We summarize the values of v^ve and 'r’eje measured at 
10 ms after the onset of merger in Table [V| This shows 
that irrespective of the models, ^ 0.3c and their 
dependence on the misalignment angle and the EOS is 
weak. Although the magnitude of Cave is quite universal, 
Ceje varies from model to model. Ceje becomes large for 
the case that the ejecta mass is large O.OSMq). This 
is consistent with the result in Fig. that the mass 

ejection proceeds in an anisotropic manner for the case 
that Meje is large. On the other hand, Ceje becomes small 
for the ejecta with small mass (< O.OIMq). This is also 
consistent with a quasi-axisymmetric ejection shown in 
Fig-IU 


comes approximately parallel with the z-axis after the 
merger. This is because the direction of the total angular 
momentum approximately preserves its initial direction, 
and the BH swallows nearly entire angular momentum 
of the system. The tilt angle of Sbh in the final state 
becomes larger for larger values of itiit,o- However, it is 
always smaller than 5°. 

The BH spin increases because the BH swallows the 
NS matter of positive angular momentum. From the 
increments of the BH mass AMbh •= ^bh f — ^bh i 


and the BH spin AaSbh I^Bny ~ ^bh,] 


we calcu¬ 


late the mean value of the specific angular momentum 
A/ := A5 'bh/AMbh that the BH gained due to the 
falling material. Here, Mbh,! is the initial BH mass and 
Sbh, pm is the value of the BH spin evaluated just before 
the merger, t ^ tmerge —1 nis. For the aligned-spin model, 
the value of Al is expected to reflect the specific angular 
momentum at the ISCO. Thus, we compare Al with the 
specific angular momentum of the ISCO of the BH with 
an effective spin parameter which we introduced in the 
previous section [see Eq. ([5^]. 

In Fig. 15, we plot A//Mbh as a function of Xeff- We 
also plot the specific orbital angular momentum at the 
ISCO of the aligned-spin BH as a function of Xeff- The 
values of A//Mbh are approximately the same as the spe¬ 
cific angular momentum at an effective ISCO. This result 
implies that the separation at which the orbital motion of 
the binary becomes unstable and the NS falls into the BH 
is given effectively by the ISCO in the equatorial motion 
around the BH with Xeff- Taking a closer look, A//Mbh 
tends to be smaller than the value for the effective ISCO. 
This is likely to stem from the gravitational-wave emis¬ 
sion that dissipates the orbital angular momentum while 
the matter falls into the BH, or the redistribution of the 
specific angular momentum due to the tidal torque. 


E. The gravitational waveform 


D. The properties of remnant BH 

We list the values of Mirr, Mbh, X, S'bh, and the tilt 
angle of the BH spin, cos~^S bu)^ at 10 ms after 
the onset of merger in Table |Vl| The mass of the remnant 
BH increases monotonically as the value of itiit,o increases 
and the EOS becomes stiff. The reason for this is that 
the tidal disruption is suppressed for larger values of itiit,o 
and stiffer EOSs, and more matter of the NS is swallowed 
by the BH. 

The dimensionless spin parameter decreases with the 
increase of itiit ,0 5 and it depends weakly on the EOS. 
In particular, for the models with qiit,o ~ 90°, the fi¬ 
nal value of the dimensionless spin parameter becomes 
smaller than the intial value of 0.75. This is because the 
dimensionless spin parameter is defined by S'bh/TTbh^ 
and the increase of is larger than S'bh for these 

models. In fact, we can see the increase of 5 'bh- The 
value of cos“^( 5 ' 3 jj/ 5 'bh) shows that the BH spin be- 


The misalignment between the orbital angular momen¬ 
tum and the BH spin causes the precession of the orbit 
and induces the modulation in gravitational waves. Also, 
the misalignment angle of the BH spin and the EOS of 
the NS affect the tidal-disruption process, and as a re¬ 
sult, gravitational waveforms are modified by them. In 
SACRA, we extract the out-going component of the com¬ 
plex Weyl scalar T 4 at finite radii and project it onto 
the spin-weighted spherical harmonic functions. Here, we 
took the axis of the spherical harmonics to be z-axis: the 
initial direction of the total angular momentum. Then 
to obtain gravitational waveforms, we integrate T 4 twice 
in time as 

t pt' 

dt' / dt"T 4 (t"). (63) 

Jo 

In the following, we plot the normalized amplitude 
Dh/rriQ or the amplitude observed at a hypothetical dis- 


h{t) = h+{t) - ihx{t) = [ 
Jo 
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FIG. 14. The volume-rendered density map of the ejecta at the time that the ejecta expands to ~ 1500 km for models H4i30 
(left), H4i60 (middle), and ALF2i60 (right), respectively. 



TABLE VI. The quantities of the remnant BH evaluated at ~ 10 ms after the onset of merger. The irreducible mass of the BH 
(Mirr,f), the mass of the BH (MBH,f), the dimensionless spin parameter (xf), the BH spin (*S'BH,f), the tilt angle of the BH spin 
(cos“^(S'gH/^BH)), and the dominant QNM frequency derived by the fitting formula Eq. ( [67| using the result for Mbh and y 
(/qnm), respectively. The results for the aligned-spin case are taken from [T7]. Eor the model APR4i90, we failed to find the 
location of the AH because of inappropriate setting of the AMR domains. 


Model 

Mirr,f[M0] MBH,f[M0] 

Xf 

8'BH,f [GMIj/c] cos 


/qnm [kHz] 

APR4iO 

6.83 

7.82 

0.85 

52 

0° 

2.56 

APR4i30 

6.93 

7.85 

0.83 

51 

< 1° 

2.49 

APR4i60 

7.16 

7.90 

0.77 

48 

1° 

2.32 

APR4i90 

- 

- 

- 

- 

- 

- 

ALF2iO 

6.78 

7.68 

0.83 

49 

0° 

2.55 

ALF2i30 

6.87 

7.76 

0.82 

50 

< 1° 

2.50 

ALF2i60 

7.14 

7.89 

0.77 

48 

1° 

2.33 

ALF2i90 

7.45 

7.95 

0.65 

41 

4° 

2.10 

H4i0 

6.74 

7.64 

0.83 

48 

0° 

2.56 

H4i30 

6.83 

7.71 

0.82 

49 

1° 

2.51 

H4i60 

7.11 

7.86 

0.77 

48 

1° 

2.34 

H4i90 

7.45 

7.95 

0.65 

41 

4° 

2.10 

MSliO 

6.74 

7.64 

0.83 

48 

0° 

2.56 

MSli30 

6.81 

7.66 

0.81 

48 

1° 

2.51 

MSli60 

7.06 

7.78 

0.76 

46 

2° 

2.35 

MS14i90 

7.44 

7.95 

0.66 

42 

3° 

2.10 


tance D = 100 Mpc as a function of approximate re- 


tarded time defined by 


Get = t- D - 2Moln^. 

Mq 

(64) 

The Fourier spectrum of the gravitational waveform 
could reflect more quantitative information. In this pa¬ 
per, we define the Fourier power spectrum of gravita¬ 
tional waves as the root mean square of two independent 
polarizations as 


(65) 

hAif) = j 

(66) 

(A = +, x). 



We will plot a dimensionless Fourier spectrum /ieff(/) •= 
fh{f) observed at a hypothetical distance D = 100 Mpc 


as a function of the frequency /, or a normalized spec¬ 
trum Dhef^{f)/mo as a function of dimensionless fre¬ 
quency ffUQ. 

Figure shows plus-mode gravitational waveforms 
of (/,m) = (2,2), (2,1), and (2,0) for models APR4iO 
and APR4i60. While the amplitude of (/,m) = (2,1) 
and (2,0) is smaller than that of (/,m) = (2,2) for the 
aligned-spin case, they could have a significant contribu¬ 
tion to gravitational waveforms for the misaligned-spin 
case. This is because the direction of the orbital an¬ 
gular momentum does not always agree with the axis 
of the spin-weighted spherical harmonic function for the 
misaligned-spin case. 

As we show in Sec. [lyg the angular velocity of the 
orbital precession is always smaller than the orbital an¬ 
gular velocity by an order of magnitude. Thus, gravita¬ 
tional waves for the misaligned-spin case have the feature 
similar to gravitational waves from the aligned-spin case 
observed from an inclined direction with respective to L 
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FIG. 15. The average specific angular momentum of material 
brought into BH normalized by Mbh, A//Mbh, as a function 
of an effective dimensionless spin parameter Xeff = xcositiit,o- 
The specihc angular momentum at the ISCO of the BH is also 
plotted. 


for each instant. Indeed, it has already shown for pro¬ 
cessing binary BH cases (see, e.g., [84]) that the wave¬ 
forms take a far simpler form in the quadrupole align¬ 
ment (QA) frame: the frame in which ^-axis agrees with 
the instantaneous direction of L. If we project gravita¬ 
tional waves onto the spherical-harmonic function in the 
QA frame, and describe these expansion coefficients by 
(/', m'), (/', m') = (2, ±2) modes are the dominant modes. 
Under the rotational transformation, these components 
mixes not only into (/,m) = (2, ±2) modes but also into 
different m modes with 1 = 2. This is consistent with 
the fact that the dominant frequency of (/,m) = (2,1) 
and (2, 0) modes agrees with the frequency of (2, 2) mode 
rather than the half (this fact is also pointed out in [36] k 
We note that different I modes do not mix under the 
rotation of the axis of spherical harmonics. 

Other m' modes, such as = (2,±1), (3,±1), 

and (3, ±3) modes, also contribute to the gravitational 
waveform. Because the phase of the m' mode is m'4>, 
where 4> is the phase of the orbit, the mixing among 
different m' modes causes modulation in the amplitude 
of the waveforms. For example, when the = (2,1) 

mode is coupled with the (2, ±2) modes, the amplitude 
exhibits modulation with the periods of 27r and Gtt in 
terms of the orbital phase, 4>. Indeed, Fig. shows 
that the amplitude observed along the 2 :-axis modulates 
primarily with the period of ~ 27r in terms of 4>. 

Obviously, the mixing of several m' modes in gravi¬ 
tational waves can occur for the aligned-spin case if we 
choose the axis of spherical harmonics which disagrees 
with the orbital angular momentum. One thing to be 
noted is that because the orbital angular momentum pre- 
cesses for the misaligned-spin case, we cannot avoid the 
situation that the orbital angular momentum disagrees 
with the axis of spherical harmonics. This implies that 
the mixing among several m' modes is unavoidable for 
the misaligned-spin case. 


Figure shows plus-mode gravitational waveforms 
(the left panels) and gravitational-wave spectra (the right 
panels) observed from different inclinations, 0 = 0°, 45°, 
and 90° with respect to the 2 :-axis. As we show in 
Sec. WA\ J is always aligned approximately with the 
z-axis, and hence 0 is regarded approximately as the an¬ 
gle between J and the direction of the observer. We take 
into account the contributions from / = 2 — 4 modes in 
plotting the waveforms. 

For MSli30, the waveforms observed from 0 = 0° has a 
similar feature to those for the aligned-spin case; the am¬ 
plitude and the frequency approximately monotonically 
increase as the system evolves. Only slight modulation 
in the amplitude and the frequency is found in the wave¬ 
forms. There are much more appreciable modulation 
in the waveforms for 0 = 45° and 90°. The amplitude 
of the waveforms becomes smaller for larger inclination 
angles. This dependence stems from the fact that for 
{V,m') = (2, ±2) modes the amplitude of the plus-mode 
waveform is proportional to (1 + cos^6>)/2. The phases 
of the waveforms agree among different values of 0. 

For MS 1190, the modulation in the amplitude and the 
frequency are larger than that for MS 1130. There is no 
monotonic dependence of the amplitude of the waveforms 
on 0. While the amplitude of the waveform for 0 = 0° 
increases as the time evolves, that of ^ = 90° approxi¬ 
mately decreases until t^et ~ 20 ms. The phase evolution 
of the waveforms is also different among different values 
of 0. This is in particular significant for the last ^ 10 ms 
of the inspiral phase. These features of the waveforms 
are due to the precession of the orbital plane. Although 
the orbital plane precesses, the angle between the z-axis 
and L is always approximately constant until the onset 
of merger for simulations in this study. Therefore, the 
angle between the direction of the observer and L is ap¬ 
proximately unchanged for 0 = 0°, and the effect of the 
orbital precession to the waveform is small. On the other 
hand, for the observer at ^ = 90°, the orbital precession 
changes the angle between L and the direction of the ob¬ 
server, and affect the waveforms. In particular, the am¬ 
plitude of the waveform is strongly suppressed when the 
orbital plane is edge-on to the direction of the observer. 

In the top-right and the bottom-right panels in Fig.[T8l 
we find that the amplitude of the spectra for MSli30 
decreases approximately monotonically for larger values 
of 0, while there is no simple dependence of the amplitude 
on 0 for MSli90. While the spectrum for the aligned-spin 
model has a flat shape (see Fig. 19 or [a [7] for the spectra 
for the aligned-spin case), there are some bumps in the 
spectra for the misaligned-spin models. In particular, 
bumps at ^ 1500 Hz in the spectra for 0 = 45° and 90° 
change the feature of the cutoff. These bumps are the 
consequences of the mixing among the different modes 
of gravitational waves. We find that I = 4 modes are 
negligible for the bumps in the spectra. Since m' = 2 
and 0 do not contribute to the modulation, {V ,m') = 
(2, ±1), (3, ±1), and (3, ±3) are primarily responsible for 
the bumps. 
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FIG. 16. The plus-mode gravitational waveforms for (/, m) = (2, 2), (2,1), and (2, 0) for APR4iO (the left panel) and APR4i60 
(the right panel). 



FIG. 17. The square of the gravitational-wave amplitude ob¬ 
served along the z-axis as a function of the orbital phase de- 
vided by 27r, 4>/27r, for models with the MSI EOS. 


The spectra of gravitational waves reflects the fate of 
the merger. In particular, as already clarified by 0111132], 
the location of the cutoff in the spectra has a correlation 
with the compactness of the NS. We expect that this also 
holds for the misaligned-spin case. However, the bumps 
of the spectra for the misaligned-spin models may make it 
difficult to determine the cutoff frequency, particular for 
the case that gravitational waves are observed from large 
inclination angles. Thus, in the following, we check the 
behavior of the bump in the spectra and the dependence 
of the cutoff in the spectra on the models qualitatively as 
a first step to understand the gravitational-wave spectra 
from the misaligned-spin BHNS. Here, we only consider 
the waveforms for an observer along the 2 ;-axis, the direc¬ 
tion of the initial total angular momentum, which would 
be the simplest case. 

We show the spectra of the waveforms for different val¬ 


ues of itiit,o with APR4 in the top-left panel in Fig. 19 
The bumps in the spectra become deeper as the value 
of ^tiit,o increases. This is because the amplitude of 
(/',m') = (2,±1), (3,±1), and (3, ±3) modes depend on 
the angle between L and the direction of the observer, 
and 0' becomes large for large values of itiit,o • For 
the models with APR4, the quasi-normal mode (QNM) 
is excited, and the cutoff frequency of the spectra re- 
ffects the frequency of the QNM. The frequency of the 
(/,m) = (2,2) least-damped mode of the QNM is given 
approximately by the fitting fomula (85] as 


/qnm = 


1 


27rMBH,/ 


[1.5251 - 1.1568(1 - Xf(.67) 


The values of the QNM frequency are summarized in 
Table |Vl| Figure shows that the spectra have a cutoff 
around the QNM frequency, although the bumps of the 
spectra for the misaligned-spin models make it difficult 
to determine the cutoff frequency accurately. 

Next, we compare the spectra of the waveforms with 
different EOSs. The top right, bottom left, and bottom 
right panels of Fig. [T^ show the spectra for itiit,o ~ 30°, 


60°, and 90°, respectively with four different EOSs. We 
again find the bumps in the spectra, and they have ap¬ 
proximately the same shape in / « 400-1000 Hz irrespec¬ 
tive of the EOS. This reflects the fact that the waveform 
does not depend strongly on the EOS for this late inspi¬ 
ral phase. The difference of the spectra due to the EOS 
becomes significant for / > 1000 Hz. The frequency of 
the spectral cutoff becomes lower in the order of APR4, 
ALF2, H4, and MSI. This is because for the stiffer EOS, 
tidal disruption occurs earlier and the waveform shuts 
down at lower frequency. Although the difference in the 
cutoff frequency is more appreciable for the spectra with 
^tiit,o ~ 30° than with itiit,o ~ 00°, we still find that the 
cutoff frequency should be, at least in principle, differ¬ 
ent among the models with itiit,o ~ 00°. For example. 
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FIG. 18. Plus-mode gravitational waveforms (left panels) and gravitational-wave spectra (right panels) observed from different 
inclination angles with respect to the z-axis. The top and bottom hgures are the results for models MSliSO and MSliOO, 
respectively. Gravitational waveforms are plotted as functions of retarded time, tret- The left axes in the plots of the waveforms 
denote the amplitude observed at a hypothetical distance D = 100 Mpc, and the right does the normalized amplitude Dhimo. 
The upper axis in the plots of the spectra denotes the dimensionless frequency, fmo, and the right axis denotes the normalized 
amplitude Dheff (/)/mo. The bottom axis in the plots of spectra denotes the frequency of the gravitational waveform, in Hz, 
and the left axis denotes the amplitude observed at a hypothetical distance T) = 100 Mpc. 


if we say the cutoff frequency as the frequency at which 
heff ~ 2 X 10“^^ is achieved, the difference of the cutoff 
frequency between models with APR4 and H4 is ~ 40% 
for both itiit,o ~ 30° and itiit,o ~ 60° cases, and the dif¬ 
ference among the EOSs is always larger than 10% for 
Aiit,o ~ 30° and 6% for itiit,o ~ 60°. By contrast to the 
cases with itiit,o ~ 30 and 60°, the difference in the spec¬ 
tra is small among different EOSs for itiit,o ~ 90°, and 
hence it might be difficult to distinguish the EOS from 
the spectra. The reason for this is that tidal disruption 
occurs so weakly that the difference in the EOS is not 
appreciable. The cutoff frequency at f ^ 2100 Hz for 
these models reflects the frequency of the QNM. Indeed, 
the value calculated by Eq. ( [67| ) agrees with this value. 


F. Emitted energy, linear momentum, angular 
momentum by gravitational waves 


A binary loses its orbital energy and orbital angular 
momentum by the gravitational radiation. The amount 
of energy and angular momentum emitted by gravita¬ 
tional waves depends on the binary parameter, and that 
is one of the interests in studying the binary merger. 
Gravitational waves could also carry linear momentum of 
the system. Non-zero linear momentum radiation causes 
recoil of the system. We evaluate these energy, linear 
momentum, and angular momentum emitted by gravita¬ 
tional waves, using the formula of [86] . 


In Table VIE we list the total energy, AE (and its ra¬ 
tio to the initial ADM mass), the total linear momentum 
normalized by the initial ADM mass, Pcw/Afo, and the 
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FIG. 19. Gravitational wave spectra observed along the z-axis of the simulation. The top-left plot shows the spectra for the 
models with APR4 and with different initial values of Ztiit,o- The top-right, bottom-left, and bottom-right plots compare the 
spectra among four different EOSs for itiit,o ~ 30°, 60°, and 9 0°, respectively. The vertical lines at the upper axis show the 
QNM frequency of the remnant BH for each models (see Table VI). 


TABLE VII. The list of total radiated energy, AE, its ratio to initial ADM mass AE/Mq, total radiated angular momentum 
AJ, total radiated linear momentum normalized by the initial ADM mass, Pcw/Afo, and the linear momentum of the ejecta 
normalized by the initial ADM mass, Peje/Mo, respectively. Eor the i90 models, the data for Peje,i was not output. 


Model 

AP[Moc^](AP/Mo[%]) AJ[GM'^/c] 

Pgw / Mq [km / s] 

Peje/Mo[km/s] 

APR4i30 

0.15(1.9) 

11 

8.1 X 10^ 

8.7 

APR4i60 

0.14(1.7) 

8.5 

5.7 X 10^ 

5.7 X 10“^ 

APR4i90 

0.098(1.2) 

6.2 

5.4 X 10^ 

- 

ALE2i30 

0.12(1.4) 

9.4 

7.8 X 10^ 

2.0 X 10^ 

ALE2i60 

0.11(1.4) 

7.8 

2.7 X 10^ 

2.0 X 10^ 

ALE2i90 

0.092(1.1) 

6.1 

3.3 X 10^ 

- 

H4i30 

0.093(1.2) 

8.5 

7.8 X 10^ 

3.3 X 10^ 

H4i60 

0.093(1.2) 

7.2 

6.3 X 10^ 

6.6 X 10^ 

H4i90 

0.085(1.1) 

5.8 

3.6 X 10^ 

- 

MSliSO 

0.074(0.93) 

7.6 

6.3 X 10" 

5.7 X 10^ 

MSli60 

0.075(0.93) 

6.6 

9.3 X 10^ 

2.8 X 10^ 

MSliQO 

0.073(0.90) 

5.5 

1.8 X 10^ 

- 
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angular momentum, A J, emitted by gravitational waves. 
We also list the recoil velocity caused by the mass ejec¬ 
tion, Peje/Mo- We take into account the contributions 
from / = 2 — 4 modes for the evaluation of the emitted 
quantities. The I = 2 modes contribute to AE and AJ 
by more than 83%, while I = 3 hy ^ 10%, and / = 4 by 
~ 3%. These fractions of the contribution depend only 
weakly on the EOS and diit,o- 

Table EUl shows that AE and AJ decrease monoton- 
ically with the decrease of the compactness of the NS. 
The same dependence of AE and AJ on the compact¬ 
ness of the NS is found for the results obtained by the 
non-spinning and aligned-spin BH-NS mergers 013, and 
it is due to the fact that a longer inspiral phase (i.e., 
longer gravitational-wave emission phase) is realized for 
a softer EOS. Table [Vn| also shows that, for a fixed EOS, 
AE and AJ monotonically decrease with the increase of 
diit,o- The dependence of AE on itiit,o is weaker than 
that on the EOS as far as tidal disruption is appreciable, 
while AJ depends appreciably on diit,o- 

The recoil velocity induced by the gravitational-wave 
emission, Pgw/TTo, decreases as the compactness of the 
NS becomes small. This is because the smaller compact¬ 
ness of the NS results in earlier tidal disruption during 
the insprial phase, resulting in an earlier shutdown of the 
gravitational-wave emission. The recoil velocity caused 
by the mass ejection is larger for the models with smaller 
compactness of the NS because the ejected mass becomes 
larger for these models. This opposite dependence of the 
recoil velocity on the compactness for Pqw and Peje re¬ 
verses the dominant component for the recoil [17] . While 
the recoil due to the gravitational-wave emission is dom¬ 
inant for models with a large compactness, the recoil in¬ 
duced by the mass ejection becomes dominant for mod¬ 
els with a small compactness. These two components are 
comparable for the case that Meje ^ O.OIMq. 


V. SUMMARY AND DISCUSSION 

We performed numerical-relativity simulations for the 
merger of BH-NS binaries with various BH spin misalign¬ 
ment angles, employing four models of nuclear-theory- 
based EOSs described by a piecewise polytrope. We in¬ 
vestigated the dependence of the orbital evolution in the 
late inspiral phase, tidal-disruption process of the NS, 
properties and structures of the remnant disk and ejecta, 
properties of the remnant BH, gravitational waveforms 
and their spectra on the BH spin misalignment angle and 
the EOS of the NS. 

We showed that a large BH spin misalignment angle 
suppresses the NS tidal disruption event by the reduction 
of the spin-orbit interaction. The remnant mass of the 
material outside the BH decreases as the misalignment 
angle of the BH increases. This dependence agrees with 
the previous results for misaligned-spin BH-NS merg¬ 
ers nsun]. Also we reconfirm the findings in 013 that 
the remnant mass increases as the compactness of the 


NS decreases. In our study, this dependence on the com¬ 
pactness of the NS is shown irrespective of the BH mis¬ 
alignment angle. The deviation of the result from the 
prediction of the fitting formula [80| is within 50% for 
Tr>AH ^ O.IMq and within 30% for M>ah ^ O. 2 M 0 , 
even though we employ a simple definition of the effec¬ 
tive spin parameter. This reconfirms the argument of [41] 
that the mass of the material outside the remnant BH 
can be modeled with a good accuracy by considering the 
result of the aligned-spin cases. 

Effects of the orbital precession are reflected in the 
tidal tail. The elevation angle of the tidal tail measured 
from the x^-plane is different for each part. Although it 
is pointed out in [41] that the elevation angle may prevent 
its elements to collide, the material of the tidal tail still 
collides with each other, and a weakly inclined torus is 
eventually formed at least for the model with itiit,o ^ 60°. 

Monotonic dependence of the remnant disk mass and 
the ejecta mass on the orbital misalignment angle was 
shown. Both M^isk and Meje decrease as the misalign¬ 
ment angle increases. However, we still found that if 
the compactness of the NS is moderate C = 0.160 or 
even small C = 0.140, BH-NS mergers with diit,o ^ 60° 
and itiit,o ^ 70° can produce disks larger than O.IMq. 
Such a system could be a candidate for the progenitors 
of sGRB. Meje > O.OIMq is achieved for diit,o < 86° with 
C = 0.140, itiit,o < 66° with C = 0.160, and diit,o < 30° 
with C = 0.176. We note that if the magnitude of the BH 
spin becomes large, more massive disk and ejecta would 
be produced. 

Eor the models with diit,o ~ 30°, the structure of the 
disk is similar to the disk formed for the aligned-spin 
BH-NS mergers. In particular, the rotational axis of the 
dense part (> 10^ g/cm^) of the disk is aligned approx¬ 
imately with the remnant BH spin. On the other hand, 
for the models with diit,o ~ 60°, we found that the axis 
of the disk is misaligned with the direction of the rem¬ 
nant BH spin initially with ^ 30°, although the mis¬ 
alignment angle of the dense part of the disk approaches 
zero in ^ 50-60 ms. While the dense part of the disk 
becomes aligned with the direction of the BH spin, an 
elevation angle of the tidal tail at large orbital radii is 
« 15° and « 30° for models with diit,o ~ 30° and 60°, 
respectively. This reflects the orbital elevation during 
the inspiral phase. It is pointed out in [82] that the 
misalignment of the disk may affect the light curve of 
the sGRB. However, since the high-density part of the 
disk with p > 10^^ g/cm^ would play main roles, the 
effect of the BH spin misalignment may not be observ¬ 
able in the sGRBs, because the dense part of the disk 
becomes aligned with the BH spin in a relatively short 
time scale.We suspect that Bardeen-Petterson-like effect 
induced by a purely hydro dynamical mechanism, such as 
angular momentum redistribution due to a shock wave 
excited in a non-axisymmetric manner of the disk, should 
work in the disk for the alignment. 

We found that the accretion time scale of the matter 
in the disk to the BH is typically ~ 100 ms, and depends 
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weakly on the binary parameters. The main mechanism 
of the accretion in the present context is the redistribu¬ 
tion of the angular momentum due to the torque exerted 
by non-axisymmetric structure of the disk, which is seen 
in Fig.|^ In reality, the viscosity induced by the magne- 
torotational instability turbulence could play an impor¬ 
tant role for this phase [88]. Since we did not take those 
effects into account, the accretion rate for the late phase 
might not be very quantitative. However, the present re¬ 
sult shows that the purely hydro dynamical effect is im¬ 
portant for the accretion of the matter, and this effect 
should be considered whenever we study the evolution of 
the accretion disk formed by BH-NS. 

We found that the velocity of the ejecta is typically 
0.2-0.3c, and has only weak dependence on the misalign¬ 
ment angle of the BH spin and the EOS of the NS. We 
also found that the morphology of the ejecta changes de¬ 
pending on the ejecta mass: Crescent-like-shaped ejecta 
with its opening angle ~ 180° is formed for relatively 
massive ejecta (> O.OSMq), while spiral-shaped ejecta 
with its opening angle larger than 360° is formed for rel¬ 
atively less massive ejecta (< O.OIMq). In particular, 
the spiral shape reflects the orbital precession. This de¬ 
pendence of the ejecta morphology was also found for 
the aligned-spin case m, and might be explained by the 
periastron advance in general relativity. 

We found that the dimensionless spin parameter of the 
remnant BH depends only weakly on the EOS of the NS, 
but it depends strongly on the misalignment angle, diit,o- 
The final direction of the BH spin becomes aligned ap¬ 
proximately with the initial direction of the total angular 
momentum. We also found an approximate relation be¬ 
tween the misalignment angle and the increase of the BH 
spin, and that the specific angular momentum that the 
BH gained during the merger approximately agrees with 
a specific angular momentum at the ISCO of the BH with 
Xeff — X^OSitiitjO • 

We showed that the mixing among the components of 
spherical harmonics occurs and causes the modulation in 
gravitational waveforms for the misaligned-spin case. In 
particular, we found that the period in the modulation 
is primarily « 27r in terms of the orbital phase. We also 
studied the dependence of waveforms on the direction of 
the observer, and found that the modulation due to the 
orbital precession becomes significant for the case that 
the observer is located along the direction perpendicular 
to the total angular momentum. The bump-shape modu¬ 
lation in the power spectrum of gravitational waveforms 
is found, and the depth of the bump becomes large as 
diit,o becomes large. 

In the presence of the bumps in the spectra, the lo¬ 
cation of the cutoff frequency becomes obscured. Nev¬ 
ertheless, for the case that gravitational waves are ob¬ 
served along the axis of total angular momentum, the 
differences of the location of the cutoff in the spectra 
among the EOSs are seen for itiit,o ^ 60°, while they are 
hardly found for diit,o ~ 90°. This result shows that, in 
principle, gravitational waves from BH-NS binaries with 


(5 = 5 and y = 0.75 contain the information of the EOS 
of the NS even if the misalignment angle of the BH spin 
is large up to ~ 60°. To discuss whether we can extract 
the information of the NS EOS from the waveform by 
the observation, we need to define a quantitative indica¬ 
tor which reflects the information of the EOS, such as a 
cutoff frequency in the spectra, in an appropriate manner 
even in the presence of the orbital precession, and discuss 
the detectability considering the noise in the signal. We 
leave these tasks for our future study. 

The dependence of the energy, linear momentum, and 
angular momentum radiated by gravitational waves on 
the misalignment angle and EOS was shown. We found 
that the recoil induced by the mass ejection dominates 
the total recoil velocity for the case that the ejecta mass 
is larger than « O.OIMq, while the recoil induced by the 
gravitational radiation is dominant for the case that the 
ejecta mass is smaller. 

Einally, we list several issues to be explored in the 
future. In this paper we studied the models only with 
(5 = 5 and y = 0.75 to focus on the dependence on the 
BH spin misalignment and the EOS of the NS. As it is 
known that the mass ratio and the BH spin magnitude 
influences on the merger process, we also need to clarify 
the dependence on these parameters with the spin mis¬ 
alignment systematically. In particular, the larger BH 
spin enhances the tidal disruption of the NS, and thus 
characteristic features of misaligned-spin BH-NS mergers 
could be revealed more clearly. Also, we plan to perform 
more detailed analysis of gravitational waveforms for the 
misaligned-spin cases, because the waveforms may con¬ 
tain rich information on the misaligned BH-NS system. 
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Appendix A: Convergence with respect to the grid 
resolution 


Table VHI compares M>ah and Meje among different 
grid resolutions for selected models. If we assume the 
first-order convergence between Y = 48 and Y = 60, 
the errors with Y = 60 results are always smaller than 
« 40% and ^ 32% for M>ah and Meje, respectively. 
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Model 

N 

M^ah[Mq] 

Meje[M0] 

ALF2i30 

40 

0.172 

3.43 X 10“^ 


48 

0.166 

3.30 X 10“2 


60 

0.156 

3.32 X 10“^ 

H4i30 

40 

0.264 

4.10 X 10“^ 


48 

0.257 

4.23 X 10“2 


60 

0.248 

4.16 X 10“2 

H4i60 

40 

0.103 

1.42 X 10“^ 


48 

9.19 X 10"^ 

1.24 X 10“^ 


60 

8.45 X 10"^ 

1.25 X 10“^ 

MSli90 

40 

2.64 X 10“^ 

8.68 X 10“^ 


48 

2.38 X 10“^ 

9.91 X lO"'^ 


60 

2.23 X 10“^ 

9.78 X lO"'^ 


Errors become large for a smaller mass. In particular, 
the error for M>ah is ~ 16% for model H4i30, while the 
error is ^ 28% for model MSliQO. 


TABLE VIIL M>ah and Meje for runs with different grid 
resolutions for selected models. 
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